Titration and hysteresis in epigenetic chromatin silencing 



Adel Dayarian 1 and Anirvan M. Sengupta 2,3 

1 Kavli Institute for Theoretical Physics, University of California, Santa Barbara, CA 
2 Department of Physics and Astronomy, Rutgers University, Piscataway, New Jersey 
3 BioMaPS Institute for Quantitative Biology, Rutgers University, Piscataway, New Jersey 

Epigenetic mechanisms of silencing via heritable chromatin modications play a major role in gene 
regulation and cell fate specification. We consider a model of epigenetic chromatin silencing in 
budding yeast and study the bifurcation diagram and characterize the bistable and the monostable 
regimes. The main focus of this paper is to examine how the perturbations altering the activity of 
histone modifying enzymes affect the epigenetic states. We analyze the implications of having the 
total number of silencing proteins given by the sum of proteins bound to the nucleosomes and the 
ones available in the ambient to be constant. This constraint couples different regions of chromatin 
through the shared reservoir of ambient silencing proteins. We show that the response of the system 
to perturbations depends dramatically on the titration effect caused by the above constraint. In 
particular, for a certain range of overall abundance of silencing proteins, the hysteresis loop changes 
qualitatively with certain jump replaced by continuous merger of different states. In addition, we find 
a nonmonotonic dependence of gene expression on the rate of histone deacetylation activity of Sir2. 
We discuss how these qualitative predictions of our model could be compared with experimental 
studies of the yeast system under anti-silencing drugs. 



I. INTRODUCTION 

One of interesting biological phenomena is the possi- 
bility for the cells with the same DNA to have different 
heritable phenotypes. Such heritable locking of cells into 
different fates without irreversible change in genetic in- 
formation is called epigenetic phenomenon pQ. One of 
the different mechanisms that can lead to epigenetic ef- 
fects is transcriptional silencing through chromatin mod- 
ification. Eukaryotic chromosomes are divided into eu- 
chromatin and heterochromatin region, based on the de- 
gree of condensation. Euchromatin regions are lightly 
condensed and genes are accessible to transcription, fn 
contrast, in heterochromatin regions, the chromosome is 
condensed throughout the mitotic cell cycle and genes are 
not normally transcribed. Consequently, the formation 
of heterochromatin is a way of silencing the expression of 
a number of adjacent genes and stabilizing gene expres- 
sion patterns in specialized cells [2] . In order for the cell 
type to be preserved in cell division, the pattern of hete- 
rochromatin and euchromatin regions has to be inherited 

The first indication for the existence of systematically 
silenced regions which are inheritable during cell division 
came from the phenomenon of position effect variegation. 
Other examples of epigenetic silencing is the HML and 
HMR Loci in budding yeast [4], and silencing of Hox 
genes, important in development of body plans, by the 
Polycomb proteins [S]. 

Different models have been proposed for silencing in 
different organisms and even for different regions of the 
genome in one organism [2 . However, there is some sim- 
ilarity among some of the proposed mechanisms [pj. In 
general, whether a region of chromosome is in the het- 
erochromatin or euchromatin state depends on the type 
of modification of histone proteins in the nucleosomes of 
the corresponding region. Here, we will discuss one of 



the models which applies to HML, HMR and telomeric 
silencing in budding yeast [4]. 

In this model, silencing initiates from a nucleation cen- 
ter which recruits certain proteins including histone mod- 
ifying enzymes. In the next step, modication of some 
of the histone tails of neighboring nucleosomes provides 
binding sites for the components of silencing complex, 
which, in turn, modify their neighboring nucleosomes. In 
such a manner, the silencing region propagates [6 a . More 
recent experiments have shown that this picture of lin- 
ear spreading of silencing from a nucleation center might 
be only valid for certain loci [7]. However, the questions 
that concern us are relevant even if there is only one re- 
gion where such spreading happens. The key question is 
what stops the spreading of the silenced region? There 
are two possible scenarios, suggested by observation from 
different loci. In some cases, there are explicit bound- 
ary elements (e.g. strong gene promoters) stopping the 
propagation [51 On the other hand, in some other 
silent regions, experiments perturbing the system by al- 
tering the abundance of pro- and anti-silencing factors 
lead to graded changes in the extent of silencing domain 
Uni E]. These observation are consistent with the 
other possibility that a self-organized stationary state be- 
tween silenced and active regions is reached, most likely 
because of the limited supply of the silencing proteins 
[T2"] . We will explore this last possibility in further detail 
here. 

Biological models of epigenetic silencing suggest 
bistable dynamics involving positive feedback loops 
where recruitment of new silencing factors is enhanced 
by the presence of chromatin bound ones in the neigh- 
borhood. Building upon that suggestion, there has been 
much computational studies of stochastic models of si- 
lencing and mean field formulation describing epigenetic 
states [I2TH7] . Earlier studies suggest that titration of si- 
lencing proteins, caused by the limited supply, has a sig- 
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nificant impact on the behavior of the epigenetic silencing 
system [12]. The purpose of this paper is to systemati- 
cally study the different regimes in which such titration 
causes qualitatively different phenomena, and explore the 
consequences and predictions of the model. 

One characteristic of nonlinear bistable systems is the 
hysteresis effects. For example, in a study of the genetic 
switch in the lac operon system, the two dimensional bi- 
furcation diagram and the corresponding hysteresis effect 
was explored by changing the abundance of two different 
molecules which affect the parameters of the system [IB] . 
In a similar spirit, it is possible to study hysteresis in epi- 
genetic silencing by exposing the cells to varying amount 
of drugs that affect histone modifying enzymes. Such 
drugs, specially histone deacetylase (HDAC) inhibitors, 
are already in use as anticancer agents [19 . In contrast to 
genetic switches like the lac system [IS], epigenetic chro- 
matin silencing has the additional feature of spreading 
along chromatin and titrating out silencing factors. This 
titration effect acts as a negative feedback competing 
with the positive feedback which gave rise to the bistable 
behavior in the first place. Therefore, it is important to 
analyze the interplay between these two phenomena. As 
we will see, depending upon the abundance of silencing 
proteins and the size of the regions affected by silenc- 
ing, one might get very different outcomes in a hysteresis 
experiment in the silencing system. 



II. MATERIALS AND METHODS 

Much is known about the silencing proteins in bud- 
ding yeast and the biochemistry of cooperative interac- 
tions between them. The key players are three proteins: 
Sir2p, Sir3p and Sir4p. These proteins form a complex 
named SIR (Silenced Information Regulator) complex 
[2D] . Sir2p is a histone deacetylase which modifies the 
neighboring histones and provides binding sites for the 
other proteins in SIR complex, namely, Sir3p and Sir4p 
[2"01 - f2"2"] . There are some other proteins which work in 
an opposing way to the silencing propagation. Partic- 
ularly, Sas2, a histone acetyltransferase, attaches acetyl 
groups to certain lysines in histone tails and prevents SIR 
complex binding [TTJ [S3] . Figure [I] shows an schematic 
presentation of the above model. 

In 1989, L. Pillus and J. Rine [1J found that in sirl 
mutants (where the nucleation effect is defective, if not 
absent), a population of yeast cells is divided into two 
distinguishable groups. In one group, HML locus is si- 
lenced similar to normal cells. In the other group, those 
loci are active and cells would not mate like a normal 
haploid. Both of the epigenetic states (silenced vs ac- 
tive) arc quite stable and are inherited most of the time 
during cell division. This observation suggests that the 
system can be thought of as being in a bistable regime, 
where two stable states can exist under the same external 
condition. In the model studied below, this bistability is 
due to competition between opposing forces and coop- 
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FIG. 1. A model of epigenetic silencing in budding yeast, 
5*. cerevisiae. In this model, bound SIR proteins helps the 
deacetylation of neighboring histone tails. Recruitment of 
additional SIR proteins to deacetylated histone tails leads to 
spreading of silenced region. 



eratively in the binding of SIR proteins. However, the 
conclusions drawn later in this study are not affected by 
the explicit mechanism of cooperatively. 



Stochastic equations describing silencing 
mechanism 



One can think of each chromosome as a 1 dimensional 
lattice, where each site corresponds to a nucleosome. 
Each site i can be in one of four possible states: 

• Bound by silencing proteins with probability Si 

• Not bound by any proteins with probability E t 

• Bound by one acetyl group with probability Ai 

• Bound by two acetyl groups with probability Di 

Figure [2] shows these possible states and the transition 
rates between them. The rate of SIR binding, which is a 
function of the ambient SIR concentration, is denoted by 
p. Free Sir2p, Sir3p and Sir4p proteins in the environ- 
ment do not form SIR complex. Instead, they form the 
complex when they are attached to a nucleosome. In the 
case where each protein is in low abundance, p is pro- 
portional to the product of the three concentrations for 
Sir2p, Sir3p and Sir4p. For our analysis, we will not need 
to know the exact form of dependence of p. We will just 
keep it as an effective parameter, monotonically increase 
with the concentration of each of SIR proteins. 

The histone acetylation rate, caused by Sas2 activity, 
is represented by a. The rate at which SIR complex falls 
off the nucleosomes is shown by 77. Also, the basal rate at 
which acetyl groups are removed from the nucleosomes is 
denoted by A. The deacetylation rate increases if adja- 
cent sites are in the silenced state. This increase is given 
by the term YijSj, where is a function of \i — j\ and 
drops significantly as this separation increases. All the 
above parameters may be position and/or time depen- 
dent. However, for the sake of brevity, this dependence 
is not explicitly written. 

We have included a double acetylation state, based on 
the fact that each nucleosome has two H4K16 sites. It 
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FIG. 2. Four possible states for each nucleosome. Chromo- 
some is modeled as a 1 dimensional lattice, where each site 
i corresponds to one nucleosome. Site i can be either in si- 
lenced (Si), unbound (Ei), monoacetylated (Ai) or double- 
acetylated (Di) state. The deacetylation rate depends on the 
silencing state of neighboring nucleosomes (the term YijSj). 



should be mentioned that there is no experimental ev- 
idence for cooperativity in the acetylation of these two 
sites. The reason we include a double acetylation state is 
to get enough nonlinearity in the system to achieve bista- 
bility (see below). This nonlinearity could be provided 
by other players and degrees of freedom not included in 
our model. Most of our conclusion regarding titration 
and bistability is independent of the particular mecha- 
nism by which this nonlinearity is introduced. 

Under the above assumptions, one gets the following 
chemical kinetics equations: 



dSi 

dt 
dEi 

dt 
dJU 

dt 
dDj 

dt 



pE l -riS l , (1) 
(rjSi - pEi) + ((A + TijS^Ai - 2aE t ), 
2a E l + 2(A + T lj S j )D t - (a + A + T^S^A,, 
aA i -2(\ + Y ij S j )D i . 



We will consider both the uniform solutions and non- 
uniform solutions in the continuum limit of these equa- 
tions. 



1. Uniform solutions 

We consider uniform steady state solutions for the set 
of equations [l] namely, we drop the subscript i and put 
the left hand sides equal to zero. Let us define 7 = S^T^. 
Using the above notation, the uniform solutions of the set 
of equations [T] has to satisfy: 

= pE - r/S (2) 

= (r/S - pE) + ((A + jS)A - 2aE) 

= 2a £ + 2(A + jS)D - (a + A + jS)A 

= aA - 2(A + jS)D 

For each set of parameters a, p and 7, we would like to 
be able to characterize how many solutions exist. The 
derivation of the solutions for the above equations is pre- 
sented in[XJ It turns out, depending on the value of the 
parameters, there can be one or two stable solutions. The 



bifurcation diagram helps to visualize different parame- 
ter values which give rise to these two different regimes 
of monostability and bistability. 



2. Bifurcation diagram 

As explained in [A| by changing the parameters con- 
tinuously, one can switch between the two regimes of 
monostability and bistability. At the point where this 
transition happens, the following equation has to be sat- 
isfied: 



7 = [S - 2S 2 



/4(1 - S)S 3 



(3) 



(2(1 - s) - v/4p-i(i-s)s) (y P s(i - s) - s) 



S(l-2S- y /4p- 1 (l-S)S^ 



(4) 



The variable S can be replaced by any value between 
and 1, as long as parameters remain positive real num- 
bers. These conditions can also be rewritten in the fol- 
lowing format: 



(a - 2(1 + aS)) , 
7= 25 ± 



(a - 2(1 + aS)) 



25" 



(1 



S 2 



(5) 
(6) 



_ 4(1 - S)S 3 
P ~ (S -2S 2 -7- 1 ) 2 ' 

The derivation of the above conditions are presented in 
|B) Equations [3] and [4] can be used to draw a plane in the 
three dimensional p - a - 7 coordinates. In fact, equations 
[5] and [6] give exactly the same plane in the 3-dimensional 
coordinates. This plane separate the the two regimes of 
monostability and bistability. It is convenient to draw the 
intersections of this plane with, for example, the constant 
p or the constant a surface. To get the former one, we 
should keep p in equations [3] and [4] constant. Instead, 
for the later case, we should keep a in equations [5] and 
[^constant. In the next section, we find it convenient to 
work with equations [5]and|6] However, for now, we stick 
with equations [3] anrT^J 

Figure [3] shows the bifurcation diagram in the a - 7 
plane (constant p). If a (acetylation rate) is relatively 
low but 7 (deacetylation effect of neighboring silenced 
sites) is high, only one solution is possible. In this case, 
most sites are in the silenced state (see [b]) . The opposite 
happens for high a and low 7 values. There also exist an 
intermediate range of values for a and 7 for which the 
system is in the bistable regime. We will denote these two 
stable solutions by Sh and Si, referring to high and low 
silencing value. At the cusp is the critical point, where 
all the solutions merge together (equation [6]) . 

For the two stable solutions at each point in the 
bistable regime, we would like to know which one is more 
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FIG. 3. Bifurcation diagram in the constant p plane, (a) The 
diagram for p = 30. The magenta line is given by equations 
[3] and [4] (b) Result of stochastic simulations. The blue and 
green points show the two monostable and bistable regions, 
respectively. The system is simulated starting from two dif- 
ferent initial conditions (high and low S). Monostable and 
bistable regimes can be differentiated based on whether the 
two initial conditions lead to one (blue cross) or two (green 
circles) different final states. 



stable. Another interesting question has to do with the 
fact that we are dealing with an spatially extended sys- 
tem. One may wonder whether it is possible to have 
different parts of the system to be in different states (si- 
lenced vs active) with an stable boundary between them. 
The subject of the next section is addressing such issues. 



B. Non-uniform solutions, coexistence of domains 
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FIG. 4. Coexistence line subdividing the bistable region, p ■ 
30. 



element stopping them from invading into each other. 
An example would be the region around the boundary of 
telomeres. 

It turns out that in the bifurcation diagram in the a - 
7 plane (constant p), the condition for the possibility of 
coexistence of domains define a line where the velocity 
of the domain boundary is zero (see [C]) . We will call this 
subspace of the parameter space the coexistence line, as 
shown in figure [4] This line starts from the cusp and 
divides the bistable regime into two sections. In the lower 
part, close to the monostable silenced regime, the silenced 
state is more stable than the active one. The opposite 
happens in the upper section. In summary, in region I or 
II of figure |1J the front between two domains is unstable 
and moves in the direction of the favorite state; the zero 
velocity condition defines the coexistence line. 

In the above discussion on coexistence of different do- 
mains, we considered a continuum system. One might 
wonder how our results would change if we had, instead, 
studied a discrete lattice model. To get insight into this, 
we simulated the stochastic system. As one may have 
expected, in the discrete version, the coexistence line 
broadens into a band of propagation failure O [25] . In 
the stochastic version of the model, within this band, 
the boundary seems to fluctuate without any noticeable 
drift. In addition, even for very large values of the pa- 
rameters (a, 7 and p), the time scale of fluctuation in the 
boundary position is quite slow. One of our future plans 
is to have a theoretical estimate on the relation between 
boundary fluctuation and the parameters of the system. 



In addition to uniform solutions discussed before, we 
would like to explore the possibility of having non- 
uniform spatial solutions for the parameter sets located 
in the bistable regime. In other words, we are looking 
for a solution which starts from one of the stable solu- 
tions (e.g. active state) and ends in the other one (e.g. 
silenced state). As we mentioned in the previous section, 
heterochromatin and euchromatin domains can occupy 
close by regions along the DNA without clear boundary 



III. RESULTS AND DISCUSSION 

A. Location of the real system on the bifurcation 
diagram 

We would like to examine whether the biological model 
of stepwise spreading of silencing fits into our mathemat- 
ical description. First, we will assume all the parameters 
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are constant. In Region I of figure [4j we do not expect 
the silenced domain to spread from the nucleation cen- 
ter. Instead, this domain should be localized around the 
nucleation center. In contrast, in Region II, the silenced 
domain spreads out from the nucleation center. Although 
this behavior is similar to the stepwise spreading model, 
it requires an explicit boundary element to stop it from 
taking over the whole active domain. 

In the telomeric regions of DNA, there does not seem 
to be an explicit boundary element stopping the spread 
of silenced domain. For example, by over expressing the 
SIR proteins, the silenced domain invades into the active 
one to some extent and then stops again [2JL This im- 
plies that the boundary between the two domains is, in 
principle, dynamic. At first glance, the stable dynamic 
boundary between the domains suggests that the system 
is actually on the coexistence line in figure [4j 

Assuming the system is on the coexistence line raises 
two concerns. The first one is that being on this line re- 
quires fine tuning of the parameter. The other issue is 
that if one of the parameters changes, e.g. p increases 
because of over-expression of SIR proteins, the system 
moves away from the coexistence line. This will cause one 
domain to invade the other one. However, in reality, this 
invasion happens only to certain extent and the bound- 
ary stabilizes at a new place. So far, our mathematical 
description does not seem to capture this behavior. 

In the above discussion, we assumed that all the pa- 
rameters are constant. In particular, the available ambi- 
ent concentrations of Sir proteins, reflected in p, was held 
constant. It turns out that by relaxing this assumption, 
not only our mathematical description explains the sta- 
bility of the boundary, but also leads to some interesting 
predictions. The details are presented in the following 
part. 

B. Consequences of limited supply of Sir proteins 

We can consider the case where the total number of 
SIR complexes, which is the sum of the proteins in the 
ambient and the ones bound to the nucleosomes, is fixed. 
In other words, there is a limited supply of SIR com- 
plexes: 

p v + Si = Stot = constant. (7) 

i 

Here, v is proportional to the volume of the system. This 
equation means that whenever a SIR complex gets bound 
to the nucleosome, the ambient concentration of available 
complexes drops. Therefore, p is a self-adjusting parame- 
ter, as opposed to being constant. We will see that there 
will be two implications from this assumption: boundary 
stabilization and coupling of different silenced regions on 
the genome. Before going forward, let us look at the 
bifurcation diagram from another angle. 

As we mentioned, the bifurcation diagram is a surface 
in the three dimensional space formed by a, 7 and p axis. 
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FIG. 5. The bifurcation diagram in the p - 7 plane (constant 
a) with a = 60. The correspondence between different regions 
of this graph and Figure [4] is shown on the picture. 

So far we have chosen to look at the intersection of this 
surface with the constant p plane (formed by a - 7 axis) . 
For the following discussion, we change this choice and 
switch to constant a plane. The diagram, which can be 
sketched using equations [JJ and [8j is shown in figure [5] 



C. Boundary stabilization without requirement for 
fine-tuning 

Let us now consider a case with a single domain bound- 
ary between the silenced and the active region. Consider 
a system located in Region II of figure [5] and assume 
there exist a small silenced domain or silencing has been 
initiated from a nucleation center. Being the favorite so- 
lution in Region II, this silenced domain invades into the 
active one. However, as silencing is spreading and SIR 
complexes get bound to the chromosome, the ambient 
SIR proteins decreases, namely, p drops. This means, 
on figure [5] the system moves vertically downward and 
approaches the coexistence line. In this way, the system 
automatically goes on the coexistence line and the two 
silenced and active domains will have a stable boundary 
between them. 

The same would have happened if we had started with 
a system in Region I with some sites in the silenced do- 
main. This time, the silenced domain would shrink and 
the system moves upward in figure [5] until it reaches the 
coexistence line. In the sense of the above discussion, 
the constraint given by equation [JJ acts as a negative 
feedback on the perturbation to the system. Figure [6] 
shows a region with a single domain boundary between 
the two epigenetic states. In this case, the constraint 
given by equation [JJ can be represented by the equation 
shown in the figure. In [D] we demonstrate how one can 
determine the location of a system with parameters S to tj 
v, 7 and L (size of the region), in the bifurcation dia- 
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FIG. 6. A silencing system with a silenced and an active do- 
mains in stable coexistence. Size of the chromosomal locus 
is denoted by L, the surrounding volume by v and the total 
number of SIR proteins by Stot- The self-adjusting parame- 
ters include the fraction of the system in the silenced state, 
x, and the density of ambient SIR proteins, p. 



gram. In other words, it is shown how to calculate the 
self-adjusting parameters p and x (the fraction of the 
system in the silenced state). 



D. Self-adjusting path in the bifurcation diagram 

The result presented so far on the stabilization of do- 
main boundary is qualitatively similar to the one studied 
in [T2] . We would like to study the effect of changing the 
parameters on a bistable silenced system and the implica- 
tions of titration of SIR proteins in greater detail. Owing 
to the constraint imposed by equation [7| altering any of 
the parameters produces an adjustment to the abundance 
of the proteins and creates a non trivial trajectory in the 
parameter space. As we will see, one gets qualitatively 
two different kinds of trajectories as the system traverses 
the bistable region. Studying the behavior of such trajec- 
tories is essential to predict and understand the results 
of a hysteresis experiment in the silencing system. 

Imagine we have a knob which allows us to play with 
the value of 7. In fact, experimentally, such a knob is 
available. By changing the concentration of nicotinamide 
(NAM), an inhibitor of Sir2p, one can effectively modu- 
late 7 [37J . We would like to know how a system changes 
as one varies the value of 7. Let us first consider the 
simple case where p is constant, as opposed to being a 
self-adjusting parameter. Figure [TJa) shows the path of 
such a system which is simply a horizontal line (magenta 
line). Figure 7j[b) shows the fraction of this system in 
the high S solution. For the case shown here, since the 
path is close to the cusp point, the size of the Region I 
and II is small. However, the shape of the curve in fig- 



ure [7lb), up to a shift along the 7 axis, is independent 
of how close or far from the cusp the path crosses the 
bistable regime. As long as we are in the monostable si- 
lenced regime or Region II (above the coexistence line) 
of the bistable regime, the whole system is in the high S 
domain. As soon as we cross the coexistence line into the 
Region I and monostable active regime, the whole system 
will be in the low S domain. 

How about when there is a limited supply of p and 
the constraint given by equation [7] is in action? Figure 
Rc) and (d) show an example. In this case, the magenta 
line in figure f7^c) and the pink line in figure |7^d) are, 
respectively, the functions p(a, 7) and x(a,j) satisfying 
the equation [2j For very large values of 7, all the sites 
are in the high S solution and the system is either in the 
monostable silenced regime (not shown in the picture) or 
Region II of the bistable regime. As one decreases 7, the 
system hits the coexistence line and the fraction x drops 
to values lower than 1. Eventually, the silenced domain 
shrinks to zero and the system enters the Region I of the 
bistable regime and then the monostable active regime. 

In the case of limited supply of SIR proteins, the sce- 
nario shown in figures [7^c) is not the only possibility and 
one can be in parameter ranges where something else 
takes place. This point is illustrated in figure [8] Fig- 
ures [8](a) shows the scenario that we already discussed. 
As shown figure [8](b) , reducing 7 causes the silencing do- 
main to shrink to its minimum size. 

The other possibility is that before the size of the si- 
lencing domain shrinks to zero, bistability is lost all to- 
gether. In this case, which happens for higher abun- 
dances of SIR proteins, the self-adjusting path passes 
through the cusp point in the bifurcation diagram (fig- 
ure |8jd)). As shown in figure |8fe), at the point where 
the bistability is lost, the size of the silenced domain is 
non-zero. 

We can consider the two silencing solutions in the 
bistable regime, Si(p(a, 7), a, 7) and Sh(p(ot, 7), ex., 7), as 
well as the single solution in the monostable regime, 
S m (p(a, 7), a, 7), along the self-adjusting path. These 
functions are shown in figure [8jc) and (f ) and calculated 
in D] Note that for the case shown in figure [8](f), the two 
solutions for SIR occupancy in the bistable regime merge 
continuously at the point where the bistability is lost. If 
the overall supply of SIR proteins is even higher, we could 
get to a regime where the system remains silenced all the 
way through the range of variation of the parameter 7 
(not shown here). 



E. Coupling different regions via ambient SIR 
concentration 



We want to analyze a situation which is inspired by 
our model system, budding yeast. Each of the 16 chro- 
mosomes in a haploid yeast has 2 telomeric regions. In 
addition, there are two silenced regions, named HML 
and HMR, located on chromosome III. The HML/HMR 
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FIG. 7. The self-adjusting path in the bifurcation diagram as 7 varies, (a) The path in the case of constant p is shown by 
the magenta line, (b) Fraction of the system in the high 5* solution for (a), (c) The path in the case of limited supply of p 
(equation [7b. (d) Fraction of the system in the high 5* solution for (c). For all panels, L = 200 sites, Stot = 110, a = 60 and 
v = l. 



loci are relatively small in size (~ 10 sites). In both 
HML/ HMR loci and telomeres, silencing is initiated by 
nucleation centers. One important difference is that 
HML/ HMR loci are surrounded by boundary elements 
stopping the silencing domain from spreading [SJ [3] ■ On 
the other hand, telomeric regions have dynamic bound- 
ary between silenced and active domains. 

Our goal is to study the effect of variation in 7 on 
this system. Since HML/ HMR loci are small in size, let 
us ignore their contribution to the constraint imposed 
by equation [7J Telomeric regions have free boundary 
and from our discussion in the previous section, we know 
how to determine the self-adjusting path in the bifurca- 
tion diagram or equivalently the self-adjusting parameter 
p(a,j). This parameter is the ambient concentration of 
SIR proteins which is also available to HML/ HMR loci. 
In other words, HML/ HMR loci read out the value of 
p{pt, 7) as it changes due to variation of 7 and the re- 
sulting effect on the state of the telomeric silenced do- 
main (see figure [9]). The possible states for HML/ HMR 
loci depends on the value of p(a, n /),a and 7. In the 
bistable regime, the two possible silencing levels are de- 
noted by Si(p(a,~f), a, 7) and Sh{p(a, 7), a, 7). In the 
monostable regime, there is only one possible solution, 
S m (p((x, 7), 01, 7). These functions are calculated in [A| 



Figure [TU] shows two examples of these functions corre- 
sponding to two different parameter regimes described in 
figure [HJ 



F. Hysteresis at HML and HMR loci 

Let us start with the following initial condition. The 
system is initially on the coexistence line as in figure |8ja). 
The silenced domain at telomeric regions coexist with 
the active domain with a free boundary separating them. 
Both HML and HMR loci are in the silenced state, i.e. 
Sh(p(cx., 7), Q!) 7)- Point a in figure 10 'a) corresponds to 



the the value of Sh{p{ot, 7), a, 7) at this initial state. 

As we decrease 7, to stay on the coexistence line, 
p(a, 7) increases. The new value of Sh(p(a, 7), a, 7) can 
be obtained as shown in [d] Point b in figure [l0[a) cor- 
responds to a value of 7 for which the system is still on 
the coexistence line. Note that, as long as the system is 
on the coexistence line, the change in Sh is continuous. 

An interesting thing happens when the path of the 
system exit the coexistence line and enters Region 
I of the bistable regime. By this time, the silenc- 
ing on the telomeric regions has shrunken to zero. 
Now, the Si(p(a, 7), a, 7) solution is more favorite than 




FIG. 8. Two different possible self-adjusting path in the bifurcation diagram as 7 varies, (a-c) The case where the size of 
silencing domain shrinks to zero as 7 decreases. In this case, the path shown in panel (a) first enters the Region I of the bistable 
regime before moving on to the monostable active region, (c) The silencing solutions along the self-adjusting path. As the 
path crosses the bifurcation diagram to enter the monostable active region, one of the two stable S solutions disappears and 
only one solution remains. There is always a discontinuity between the two stable solutions in the bistable regime, (d-f) The 
case where the bistability is lost before the size of the silencing domain shrinks to zero, (f) As the path crosses the bifurcation 
diagram through the cusp point to enter the monostable active region, the two stable solutions merge in a continuos manner 
and form the single monostable solution. For all panels, L = 500, a — 60 and v = 1. In panels (a-c), Stot = 110 and in panels 
(d-f), Stot = 180. 



Test region (e.g. HMR) 



Boundary 
element 




FIG. 9. Several silencing systems coupled through a shared reservoir of ambient SIR proteins. Perturbing one of the parameters, 
e.g. 7, leads to the self-adjustment of the boundary location in regions with free boundary, such as telomeres. Consequently, the 
concentration of ambient SIR proteins, p, changes. This, in return, affects the silencing in other regions with explicit boundary 
elements such as HMR. 



Sh(p(a, 7), a, 7). Therefore, the state of HML/HMR loci 
would change to the lower, more active solution. This 
transition is shown in figure |10[ a) by the pink arrow. 
The new silencing state of the HML/HMR loci is repre- 
sented by point c in the figure. As we keep decreasing 7, 
the Si(p(a, 7), a, 7) solution decreases as well. Eventu- 
ally, the system crosses the bifurcation line (point d) and 
goes into the active monostable regime (point e). 

What happens if we increase 7? From point e to d to c 



in figure [T0| the state of HML and HMR loci goes back on 
the same path as before. However, at point c, where the 
system hits the coexistence line, the level of silencing at 
HML/HMR loci takes a new path. Previously, when we 
approached this transition point from right, HML/HMR 
loci were in the Sh{p(a, 7), a, 7) solution, whereas this 
time, they are in the Si(p(a, 7), a, 7) solution. If we in- 
crease 7, the state of these loci will stay on the lower 
branch and move towards point /. One counterintuitive 
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FIG. 10. Hysteresis effect in the silencing level at HML/ HMR. (a) Curves represent the silencing solutions along the self- 
adjusting path shown in figure |8^a). Silencing level at HML/ HMR is always located on the blue curves. As 7 is decreased, the 
silencing at HML/ HMR follows the points a — s> b — > (with a jump) c — > d — > e. Subsequent increase in 7 changes the silencing 
value ase->d->c->/. (b) similar to (a) but for the path in figure |8|d). Note that the transition between points b' and 
(c', d!) is continuous. Moreover, starting from point e', increase in 7 can take the system to either the lower (point /') or upper 
(point b') branch. 



behavior in the above discussion is that, at point / com- 
pared to point c, although 7 is higher, the silencing has 
reduced. The same was also true between points c and 
d: 

fed 

Si (p(a,7),a,7) < S t (p(a, 7), a, 7) < S l (p(a, 7), a, 7). 

(8) 

Let us now examine the parameter regime where the 
self-adjusting path is described by figure [sjd), the above 
picture would get modified. The silencing solutions (Sh, 
Si and S m ) for this case are depicted in figure |l0[ b). 
Again, assume that we start at point a' in figure \l0\h) 
where the silenced domain at telomeric regions coexist 
with the active domain with a free boundary separating 
them and both HML and HMR loci are in the silenced 
state, i.e. Sh(p(a, 7), a, 7). When we decrease 7, the si- 
lencing solution at HML/ HMR loci moves to point b' in 
figure 10 'b). As the path of the system exit the coexis- 



tence line and the bistable regime through the cusp point 
(see figure [8^d)) , the silencing solution at HML/ 'HMR 
loci changes smoothly and reaches the point denoted by 
(c', d!) in figure [Toj^ b) . Note that, in contrast to the case 
shown in figure TTOT a), there is no jump in the silencing 
level at HML/ HMR loci. Further decrease in 7, will even- 
tually move the silencing level towards the monostable 
regime (point e'). 

As we start to increase 7, the silencing level at 
HML/ HMR loci goes from point e' back to point (c', d'). 
At this time, with further increase of 7, the self-adjusting 
path in figure [8jd) crosses the cusp and moves along the 
coexistence line. As we mentioned before, the two stable 
solutions in the bistable regime are equally stable when 
the system is on the coexistence line. As a result, and in 
contrast to the case described in figure 10 'a), the silenc- 
ing level at HML/ HMR loci does not necessarily follow 
the active branch (lower branch) towards the point /' 
in figure |10[ b) . Instead, it can go back from the upper 
branch towards the point b' with equal probability. In 



other words, by sweeping 7 up and down, one does not 
necessarily obtain the familiar hysteresis loop. 

The above discussion was based on a deterministic de- 
scription of the system. However, that discussion still 
guides our conclusions for the behavior of a real system 
with stochastic effects. If we lower Sir2 activity to go 
to the monostable active region, and then slowly bring 
the activity up, at some point, we come to the bound- 
ary of the bistable region. In case the system is passing 
through the cusp, one expects that the population breaks 
into two subpopulation, one following the upper branch, 
while the other takes the alternative lower branch. As 
long as the relative weights of the two subpopulation be- 
tween these two branches remain different for different 
initial conditions, we still see a hysteresis effect. 



IV. CONCLUSION AND OUTLOOK 

In this paper, we concern ourselves with hysteresis ef- 
fect in a model of chromatin silencing in budding yeast. 
We compute the bifurcation diagram of the system and 
analyze the result of sweeping up and down one of the 
parameter, namely, the catalytic activity of Sir2. In par- 
ticular, we consider the case where there are limited sup- 
ply of SIR proteins and show that the resulting depletion 
effects could alter our usual view of hysteresis. Two most 
interesting observations are the followings. 

The first is that, as the Sir2 activity (represented by 7) 
is reduced, it is possible to go from the silenced branch 
of the system to the active branch without a jump. The 
other counterintuitive observation is that, along the ac- 
tive branch in the bistable regime, the silencing decreases 
as the Sir2 catalytic activity increases. Both of these ob- 
servations lead to qualitative predictions that could be 
compared with experimental results. Features like this 
are characteristic properties of switches based on chro- 
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matin modification where the extent of the silenced do- 
main is a self-adjusting variable and stands in contrast 
with more familiar examples of switches based on tran- 
scriptional regulatory feedback. 

One available set of tools from biomolecular chemistry 
is a large number of drugs than can alter the activity of 
histone modifying enzymes I28--31]. In fact, the catalytic 
activity of some of the relevant proteins in the yeast si- 
lencing system, such as Sir2 and Sas2, could be affected 
by such chemical drugs. Of special interest are HDAC 
inhibitors of Sir2 like nicotinamide (NAM) [57] and split- 
omycin |32j . By monitoring the presence and absence of 
long term epigenetic memory as we sweep up and down in 
the parameter space using such inhibitors, we should be 
able to directly test the predictions made in this study. 



whereas for relatively high values of a, the solution is at 
low 5 (figure [Tl| b)). There is also a regime of param- 
eters where there are three real solutions (figure [TT|c)). 
The middle solution is unstable, whereas the other two 
solutions at low and high values of 5 are stable. We will 
denote these two stable solutions by 5/ and Sh, respec- 
tively. When the parameters allow us to have two stable 
solutions, we are in the bistable regime. As we play with 
the parameters, for example by increasing/decreasing 7, 
two of the three solutions merge together (figure [TT|d)) 
and we are on the boundary between the bistable and 
monostable regime. At this point, the curves are tangent 
to each other. The critical point or the cusp in the bi- 
furcation diagram is defined as the point where all three 
solutions merge together. 



V. ACKNOWLEDGEMENTS 



Appendix B: Derivation of bifurcation diagram 



We thank Vijaylakshmi Nagaraj and Mohammad 
Sedighi for many useful discussions. AD was supported 
by NSF PHY11-25915 while AMS acknowledges support 
of National Human Genome Research Institute grant 
R01HG003470. AMS thanks the Kavli Institute for The- 
oretical Physics for supporting his visit during which part 
of this project was completed. 



Appendix A: Derivation of uniform solutions 

By eliminating the variables in the equation[2] we find: 



E=—- A = 2ar]S 

P ' P (A + 7 5) '' 



D = 



Let us also rescale the parameters by A as follow: 

1 P 1 a 1 7 

P = - " = v! 7 = -■ 
r) A T) 



(1) 



(2) 



From now on, we drop the primes on the rescaled param- 
eters. Since the sum of the probabilities has to be 1, we 
have: 



, , , 1 2 a 

5 1 + - + 



1; 



which we can rewrite as: 



5 = 



P (l + 7 5) 2 



[(1 + p) (1 + jS) 2 + 2a(l + jS) + . 



(3) 



Figure [TTJ shows the graph of left hand and right hand 
side of the above equation as a function of 5 for a few dif- 
ferent combination of the parameters. The above equa- 
tion is of degree 3 and can have at most three real roots. 
For certain sets of parameters, referred to as the monos- 
table regime, there is only one real solution (figure 11 a) 
and (b)). For relatively small values of a (or high values 
of 7 or p), this solution happens at high 5 (figure 11 'a)), 



In [A] it was shown that by changing parameters con- 
tinuously, one can switch between the two regimes where 
the set of equations [2] have one or three real solutions. At 
the transition between these two state, the two curves in 
figure [TT] touch each other at a point. This is the point 
where two of the solutions merge and disappear or a de- 
generate solution appears and eventually give rise to two 
solutions, depending on the direction that we change the 
parameters. At this point, not only the equation[3]is sat- 
isfied, but also the derivative of both sides with respect 
to 5 should be equal. Let us rewrite equation [3] as: 

p (1 + 1S) 2 = 5 [(1 + p) (1 + 7S) 2 + 2a(l + jS) + a 2 ] . 

(1) 

Putting the derivative of both sides equal, and using 
equation [JJ we get: 

2p 7 5 (1+75) = p (l+ 7 5) 2 +5 2 [2 7 (1 + p) (1 + lS) + 2a 7 ] . 



This implies for the transition point: 



a=(l + 7 5) 



p(l-S) _ p(l + 7 5) 
5 2 7 s 2 



(2) 



(3) 



After replacing a in equation [JJ by the above equation 
and dividing both sides by 5(1 + 75) 2 , we get: 



gft - S) 

5 



p(l-5) p(l + 7 5) 



5 



p(l-5) 



2 7 5 2 

p(l-S) P (l + jS) 
5 2 7 5 2 ' 



The reason we take the positive root is because a > 0; 
therefore, the term in the bracket in equation [3] is posi- 
tive. We can solve the above equation for 7: 



7 = 5 - 25 2 - 



'4(1 -5)5 3 



(4) 
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FIG. 11. The intersections of nullcline curves. Graphs shows the left hand side (magenta) and the right hand side (blue) of 
equation[3] for different sets of parameters, (a) Monostable silenced regime, a = 30, p = 30 and 7 = 40. (b) Monostable active 
regime, a — 115, p — 30 and 7 = 40. (c) Bistable regime with two stable solutions and one unstable solution, a — 85, p — 30 
and 7 = 40. (c) System on the bifurcation line between the bistable regime and the monostable active regime, a = 97, p = 30 
and 7 = 40. 



We can replace the above in equation [3] to get: 

(2(1 - 5) - y/4p-i(l-S)s) (y/pSQ. - S) - s) 

o = — — — . 

s(l-2S- ^4p- 1 (l-S)S^j 

(5) 

In equations [4] and [5j for each value of p, 5 can take 
any value between and 1 as long as both a and 7 are 
positive real numbers. 

There is one more case that we did not mention and 
is not shown in figure [TT] It is possible to have a situa- 
tion where all three solutions merge together. This case 
is similar to figure |TT|d) , with the difference that the 
two curves intersect only at one point. To be in such a 
situation, in addition to equations [T] and [2j the second 
derivative of both sides of the equation [T] with respect to 
5 has to be equal. Putting the second derivatives equal, 
and using equation [l] and equation [2j we get: 



7/ o - 2a - 2(1 + p) 
37(1 + P) 



(6) 



The subscript C is meant to indicate the critical point or 
the cusp. Note that the parameters in the above equation 
have to satisfy equations [T] and [2] as well. 

Equations Hand [5] are the consequence of the two con- 
ditions [T] and [2] that we have imposed. Instead of solving 
for a and 7, we could have chosen to solve for a and p, 



or p and 7. If we solve for p and 7 in equations [T] and [2] 
we find: 



(a - 2(1 + aS)) , 
7 = — ± 



P = 



25 

4(1 - S)S 3 
(S-2S 2 -~/- 1 f 



(a - 2(1 + aS)) 



1 2 



25 



(! + "'■ 



■7) 



(8) 



where 7 in the second equation has to be replaced from 
the first one. 



Appendix C: Derivation of non-uniform solutions 

Given that we are dealing with a spatially extended 
system, in addition to uniform solutions for equation [2] 
we would like to explore the possibility of having non- 
uniform spatial solutions for the parameter sets located 
in the bistable regime. For a system with N nucleo- 
somes, there are 4 N possible distinct states. We can- 
not directly solve the time-independent solutions of the 
stochastic system given by the set of equations [T] There- 
fore, we will resort to the continuum limit approximation. 
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In this limit, set of equations [T] become: 
dS(x 



dt 
dE(x) 
dt 



dA{x) 
dt 



P E{x) - V S(x) , (1) 
r}S(x) - pE(x) - 2aE(x) + 
(^X + J r(x-y)S(y)dy^ A(x) , 

2a E(x) + 2 ^A + J T(x - y)S(y)dy^j D(x) - 

a + X + J T(x - y)S(y)dy^j A(x) , 

aA(x) - 2 ( A + [r(x - y)S(y)dy] D{x). 



dDjx) 
dt 

We can Taylor expand S(y) around x. Since T(x — y) 
falls sharply as \x — y\ increases, we will only keep up to 
the third order in the expansion: 

r(rr - y)S{y)dy = J T{x - y) (s(x) + (y - x) ^ + 
(y - x) 2 d 2 S(x) 



jS(x) + 72 



dx 2 
d 2 S(x) 
da; 2 



dy 



(2) 



The second term in the Taylor expansion disappears since 
T(x — y) is symmetric. We have also defined: 

1=1 F(x-y)dy & j 2 = / T(x-y) ^ ^ dy (3) 



Replacing [2] in the set of equation [T] and simplifying the 
equation we arrive at: 
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d 2 S(x) 



dx 2 
If we define: 

V(S) = +s- 



= — 1 — iS(x) + a 



S(x) + y/pS(x)(l-S(x)) 



p(l - S(x)) - S(x) 



(4) 



1 S 2 



S' + y/ P S'{\ - S>) 



dS'. (5) 



2~ "7 p(i-S')-S' 

then equation [4] can be written in the following form: 



72- 



d 2 S(x) 
dx 2 



dV(S) 
dS ' 



(6) 



The similarity between [6] and the formula for the motion 
of a particle in a potential field in classical mechanics is 
clear. In this analogy, S here plays the role of the position 
of the particle and x plays the role of time. 

For the two uniform stable solutions of Equations [T] 
Si and Sh, the right hand side of equation [4] is zero. At 
those points, the potential V, defined in equation [5j is 
flat (^7 |s, = ^-\s h =0). Using equation pi we can nu- 
merically calculate the values of V for the points between 



numerical integration for different combination of param- 
eters within the bistable regime. 

We are looking for a solution which starts from one 
of the stable solutions (e.g. Si) and ends in the other 
solution (e.g. Sh)- From our experience in classical me- 
chanics with equations in the form of [6j we now that the 
necessary condition is (Figure [l2jb)): 



(7) 



V(Si) = V(S h ) 



It turns out that in the bifurcation diagram in the a - 
7 plane (constant p), this condition defines a line that 
we will call the coexistence line. Figure [4] shows this 
coexistence line. Note that we can use this potential 
only to describe the zero-velocity fronts, and the general 
traveling solution. 



Appendix D: Configuration of a system with limited 
supply of SIR proteins 

If the system is in the bistable regime, there are two 
possible states, Si(p, a, 7) and Sh(p, a, 7). In the silenced 
monostable or active monostable region, there is only one 
possible solution, S m (p,a,-f). Note that for fixed a and 
7, these solutions are monotonically increasing function 
of p. 

Consider a particular value of 7 and a. Assume this 
value is high enough so that for certain range of p the sys- 
tem is in the bistable regime. Lets consider the following 
function: 



(1) 



k = S m (p,a,j) L + p v 

if p in active monostable region, 

h = Si(p,a,j) L + p v 

if p in Region I of bistable regime, 

^>3 = Si{p, a > 7) (1 - L + s h(p, a,~f) x L + p v 

if p on the coexistence line and where < X < 1, 

k = S h (p,a,-f) L + pv 

if p in Region II of bistable regime, 

^5 = S m (p,a,-f) L + pv 

if p in silenced monostable region. 



The variable x represents the fraction of the system in 
the Sh domain. It can be different from zero or one only 
for a particular value of p for which the system is located 
on the coexistence line (See figure 13). Note that is a 
monotonically increasing function of p and x. 

For a fixed value of 7, a and Stot, to determine what 
the configuration of a system is and where in the bifur- 
cation it is located, one has to find: 



the two stable solutions. Figure 12 shows the result of 



pia, 7) and xia, 7) such that <pip,a,^,x) = Stot- (2) 
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FIG. 12. Potential V(s) for different combination of parameters in the bistable regime. At each point in the bistable regime, 
there are two stable solutions for Equation [l] Si and Sh- The graph is only drawn for values of S between these two solutions, 
(a) q = 90, p — 20 and 7 = 55. (b) a — 90, p — 20 and 7 = 50. (c) a = 90, p = 20 and 7 = 49. The correspondence with figure 
PA is as follow: panel (b) corresponds to the coexistence line; panel (a) is located in Region II; panel (c) is in Region I. 



Ambient SIR protein v Total SIR protein 



/ 

S h (p,a,-y)xL + Si(p,a,-y) (1 - x) L + pv = S u 

V 



Silenced domain 




Active domain 



Chromosome position 

FIG. 13. Constraint imposed by limited supply of SIR pro- 
teins. 



Note that, we should have really written p(a, 7, Stot, L) 
and x(a,j,Stot,L), however, for the sake of brevity, we 
did not write the last two parameters. How can we calcu- 
late p(a, 7) and x(a, 7)? Let us consider some particular 
values of p which are of interest. In the active monostable 
region, the minimum value of p is and the maximum 
value happens when we touch the bifurcation line (green 
line) in figure [5] from below. We call this value /0fe«(a,7) 
(b stands for bifurcation and u for active). Let us refer 
to the value of p on the coexistence line by p z (a,~/). As 
we increase p, we hit the green line again, this time on 
the boundary between bistable regime and the silenced 
monostable one. Let us call this value pb s (a, 7) (b stands 
for bifurcation and s for silenced). With this notation, 
we have: pb u < Pz < Pbs- Using this notation and the 
definitions given in equation [l] we have: 

fa(p b u(a,-f),a,j) < fa(p z (a,-f),a,j,x = 0) < 
fa(Pz(a,j),a,j,x = 1) < fa(p bs (a,j),a,j). (3) 

Also, note that: 

0i(p6u(a,7),a,7) = fa(pbu(a, 7), a, 7) and 
(j} 4 (pbs(a,-f),a,j) = fa(p bs (a,j),a,-f). ( 4 ) 



The first step in determining the configuration of a 
system and where in the bifurcation diagram it is lo- 
cated is to compare S to t with the 4 values involved in 
equation |3] If Stot < 4>i(pbu(ct, 7), a, 7), the system is 
located in the active monostable region. Similarly, if 
Stot > 4>b(,Pbs(p<-i 7): ol, 7), the system is located in the si- 
lenced monostable region. If (f>^pb u {a, 7), a, 7) < S to t < 
4>3{Pz{o>>i 7)) a, 7; x = 0), the system will be in Region I of 
figure[5] On the other hand, if (f>3(p z (a, r y),a,j,x= 1) < 
Stot < 05(pbs(a, 7), a, 7), the system will be in Region II. 
For each of these regions, one can numerically solve the 
corresponding fa in equation [l] for different values of p 
and find the one that satisfies the constraint given by 
equation [2j 

The only remaining case is when 

fa(Pz(a,-f),a,j,x = 0) < Stot < fa(pz(a,j),a,j,x = 1), 

which corresponds to a system with domains of silenced 
and active regions coexisting with each other. The frac- 
tion of the system in the Sh domain is determined by 
satisfying: 

Si(p e (a,'y),a,j) (1 - x) L +S h {p z (a,j),a,j) x L + 

p x {a,i) v = Stot, (5) 

which implies: 



x _ (Stot - Pz(a,l) v ~ Si(p z (a,j),a,'y) L) 
(S h (p z {a,j),a,j) L — Si(p z (a,j), a, 7) L) ' 



In summary, we showed how to calculate the self- 
adjusting parameter p and the configuration of the cor- 
responding system. 
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